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Using results from colloid science we derive interaction 
potentials for computer simulations of mixtures of soft or 
hard ellipsoids of arbitrary shape and size. Our results are 
in many respects reminicent of potentials of the Gay-Berne 
type but have a well-defined microscopic interpretation and 
no adjustable parameters. Since our potentials require the 
calculation of similar variables, the modification of existing 
simulation codes for Gay-Berne potentials is straightforward. 
The computational performance should remain unaffected. 



I. INTRODUCTION 

In molecular simulations [1,2] short-range attractive 
and repulsive interactions are typically represented us- 
ing Lennard-Jones (LJ) 6-12 potentials: 
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where a/r is the dimensionless ratio of the effective par- 
ticle diameter and the interparticle distance. While the 
r~ 6 part has a physical origin in dispersion or van der 
Waals-interactions, the r~ 12 repulsion is chosen by math- 
ematical convenience. For large molecules the evaluation 
of the interaction potential involves a computationally 
expensive double summation of Eq. (1) over the respec- 
tive (atomic) interaction sites 
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or the evaluation of a double integral 
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in the corresponding continuum approximation for bod- 
ies with simple geometric shapes and number densities 
Pi{r) of interaction sites. We will refer to interaction en- 
ergies obtained by (numerically) evaluating Eq. (3) as the 
"Hamaker" potential. 
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As an alternative Gay and Berne [3] (GB) proposed 
the use of more complicated single-site interaction po- 
tentials for rigid molecules. Their approach is based on 
a heuristic modification of a Gaussian overlap potential. 
While GB potentials provide a computationally efficient 
way to introduce anisotropic interactions in numerical 
studies of liquid crystalline systems [4,5], they have fre- 
quently been criticised for their unclear microscopic in- 
terpretation [6]. In the present paper, we use results from 
colloid science [7] to derive approximate interaction po- 
tentials for mixtures of ellipsoids of arbitrary size and 
shape which have a well-defined microscopic interpreta- 
tion and no adjustable parameters. 

The paper is organized as follows: After introducing 
the Gay-Berne potential in section II, we review in sec- 
tion III the Hamaker theory for two spheres of arbitrary 
size and develop a relatively simple approximation of the 
interaction potential which is valid at arbitrary distances. 
In section IV we generalize this expression to the case of 
interacting ellipsoids. Section III presents a numerical 
test of this approximation for the case of pole contacts 
between aligned ellipsoids. In section VI we suggest com- 
putable expressions for the orientational dependence of 
the interaction potential which are in many respects rem- 
inicent of those familiar from GB potentials. A numerical 
test of the proposed interaction potential for ellipsoidal 
particles of different shape at arbitrary relative position 
and orientation is presented in in section VII. We con- 
clude with a brief summary in section VIII. 



II. THE GAY-BERNE POTENTIAL 

A rigid body i is specified by its center position fj, its 
orientation (expressed, for example, via a rotation matrix 
Aj for the transformation from the lab frame to the body 
frame) and its shape. In the case of ellipsoids, the shape 
is given by three radii ai,bi, Ci which can be used to define 
a "structure matrix" 
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in the body frame of the ellipsoid. 

The most general form of a Gay-Berne potential for 
dissimilar biaxial ellipsoids was introduced by Berardi, 
Fava and Zannoni (BFZ) [8] as a product of three terms: 

U(Ai,A 2 ,fi 2 ) = C/ r (Ai, A 2 ,fi 2 ) 

xt ? i 2 (Ai,A 2 ) xi 2 (Ai, A 2 ,fi 2 ) (5) 
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where (ri 2 ) ^12 is the (unit) vector between the center 
positions: F12 = r 2 - r 1 and fi 2 = ri 2 /|ri2|. 

The hrst term controls the distance dependence of the 
interaction and has the form of a simple LJ potential 
Eq (1) 
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where the interparticle distance r\i is replaced by the 
distance h\i of closest approach between the two bodies: 



/112 = min(|rj - fj\) V(i e Body 1, j e Body 2) 



(7) 



The position of the potential minimum, (2 1 / 6 — 7)0" is 
shifted empirically relative to the Lcnnard- Jones value 
2 1 / 6 cr. Typically 7 = 1. The well depth is mm(U T ) = 
£gb- 

In general, the calculation of /i 12 is non-trivial. For 
ellipsoids a suitable scheme was worked out by Perram 
et al. [6,9,10]. These authors also clarified the meaning 
of the distance 
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which is usually employed together with the Gay-Berne 
potential [11]. Eq. (8) is an approximation which fails, 
for example, in the case o f two spher es with unequal radii 
a ± < a 2 where 0^2 = \/2(a\ + of) w \/2a 2 > ai + a 2 . 
In this article we always use the correct contact distance 
hi2- Fig. 4a and b provide a comparison of the quality 
of the various approximation schemes, if h\2 is replaced 
by the Gay-Berne approximation hf 2 B . 

The two other terms in Eq. (5) control the interaction 
strength as a function of the relative orientation and po- 
sition of the ellipsoids. The second term [12] introduces 
an empirical exponent v. 
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The third term has the form 



Xl2 (A 1 ,A 2 ,f 12 ) = [2ff 2 Br 2 1 (A 1 ,A 2 ) fi 2 ]". 
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with 



and 



Bi 2 (Ai, A 2 ) = AfEiAi + A 2 r E 2 A 2 



E, 




(15) 



where e tt j, e^, e c % characterize the relative well depth for 
side-to-side, face-to-face and end-to-end interactions be- 
tween two ellipsoids of type i. fi is another empirical 
exponent. 

To summarize, the physical problem of a mixture of 
colloidal particles of equal composition, but of different 
sizes and (ellipsoidal) shapes, is defined via Eqs. (1) and 
(3). It requires the specification of a material constant 
and of the shapes ai,bi, Ci of the involved particles. Gay- 
Berne potentials introduce additional adjustable param- 
eters: the shift parameter 7, the empirical exponents v 
and /1, and three energy parameters per particle type. 
These parameters are usually adjusted by fitting Eq. (5) 
to the numerical evaluation of Eq. (2) for small assemblies 
of suitably arranged Lennard- Jones particles [3] or spe- 
cific organic molecules [13]. Note, that there are no ad- 
ditional parameters specifying the interactions between 
ellipsoids of different shape. Rather, Eqs. (6) to (15) 
provide heuristic "mixing rules" for this case. 

In the following we will partially justify the orientation 
dependent part of the Gay-Berne potential and the im- 
plicit mixing rules for particular choices of the adjustable 
parameters. Note, however, that the product ansatz of 
Eqs. (5) incorrectly reduces the attractive and repulsive 
parts of the interaction between extended objects at ar- 
bitrary distances to simple power laws with distance in- 
dependent, shape and orientation dependent prcfactors. 
To overcome this problem we will abandon the strategy 
initiated by Gay and Berne who sought modifications to 
the Lennard-Jones potential for point particles. Instead 
we will try to preserve the case of interacting spheres of 
finite volume as a proper limit. 



III. HAMAKER THEORY 



Eq. (3) can be solved exactly for two spheres of radius 
»i < 0-2, volume Vi = -j-af at a distance ri 2 = (a\ + 
0-2) + h\2 with h\2 > 0. For the attractive part of the 
interaction Hamaker [14] obtained 
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where ^4i 2 is usually referred to as Hamaker's constant. 
Using LJ units ^4i 2 is given by Ai 2 = 4ir 2 e(pa 3 ) 2 . Simi- 
larly, we found for the repulsive part of the LJ potential: 
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FIG. 1. Attractive part Ua of the potential energie multi- 
plied by the inverse asymptotic distance dependence h 6 as a 
function of the distance of closest approach h for pole con- 
tacts between differently shaped ellipsoids: (a) Hamaker po- 
tential, (b) the approximation proposed in the present paper, 
and (c) a Gay-Berne 6-12 potential adjusted to reproduce 
the energy minima within the Deryaguin aproximation. The 
curves converge at large distances, because we have used pro- 
late ((a, b, c) = (1,6, 6)cr), oblate ((a,b,c) = (2,2, 9)<j) and 
spherical (a — 36 3 <r) ellipsoids of identical volume. The 
plots contain results for sphere-sphere (thick solid line), pro- 
late-prolate (gray lines), oblate-oblate (narrow dark lines), 
and prolate- oblate (dotted lines) contacts. 
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Some insight can be gained by considering three limit- 
ing cases: (i) distances which are smaller than the (cur- 
vature) radii of the spheres, (ii) a small sphere (i.e. a 
point particle) at an intermediate distance from a much 
larger sphere, and (iii) the large distance limit of a 6-12 
Lennard-Jones potential with appropriately renormalized 
prefactor: 
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Fig. la shows a log-log plot of UaJt- 6 which illustrates the 
deviations from the asymptotic power law at small dis- 
tances. The figure also contains numerical result for the 
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FIG. 2. Sphere size dependence of the position (a) and 
depth (b) of the energy minimum: Hamaker potential (•) 

Eqs. (16) and (17), our approximation ( ) Eqs. (20) and 

(21), large/Deryaguin ( ) and small/Lennard- Jones (■ ■ ■) 

sphere limits Eqs. (18) and (19). 



attractive part of the Hamaker potential for pole contacts 
between prolate and oblate ellipsoids. Qualitatively, the 
curves resemble each other. They converge at large dis- 
tances, because the particles were chosen to have iden- 
tical volumes. In contrast, at small distances the inter- 
action strongly depends on the relative orientation and 
position of the non-spherical particles. 

Since Eqs. (16) and (17) are too complicated for an 
approximate generalization, we have instead developed a 
suitable combination of the three limiting cases discussed 
above: 



36 \ 01 + a 2 n 12 

3 



«i 



a 1 



Ur 



AV2_ 

2025 



h) 



W2 

6 

' 1 



"2 



a-2 + /112/2 
45 2ai a 2 1 

0-2 hi2 



a 1 



56 a\ 

3 



a x +/ii 2 /60 1 /3; \a 2 + h 12 /60 1 / 3 



"2 



(20) 



(21) 



Prefactors were chosen in such a way that Eqs. (20) 
and (21) reproduce the limits Eqs. (18) and (19) of the 
exact results. As a consequence, the approximation is 
fairly reliable on all length scales (compare, for exam- 
ple, the thick solid lines in Figs, la and b). This is also 
demonstrated in Fig. 2 where we show the sphere size 
dependence of the depth and position of the minimum of 
Ur + Ua- While the short distance expansion becomes 
reliable for sphere radii a > 5<t, Eqs. (20) and (21) es- 
sentially reproduce the exact results for arbitrary sphere 
sizes. 

As pointed out above, potentials of the Gay-Berne type 
Eq. (5) cannot describe the complex distance dependence 
of the interaction. But what about the potential mini- 
mum? Within the short distance/large sphere expansion, 
its depth U m i n and position h±2, m in are given by: 
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For the Gay-Berne potential, on the other hand, 
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in the distances in Eq. (6) is thus a natural consequence 
of the insistence on a 6-12 potential. Furthermore, the 
comparison suggests a relation between the energy scale 
of the GB potential and the Hamakcr constant 
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as well as the choices /x = 1 and e, = a^/cr. Note, that 
this result represents a first justification for one of the 
empirical GB mixing rules. 



IV. INTERPOLATING BETWEEN THE 
DERYAGUIN APPROXIMATION AND THE 
LARGE DISTANCE LIMIT 

In the general case of ellipsoids of arbitrary shape, rela- 
tive position and orientation, the small and the large dis- 
tance limit remain (qualitatively) unchanged. However, 
for strongly non-spherical ellipsoids with a, <C there 
are new regimes for intermediate distance Oj -C h -C c,. 
For example, one finds for two thin prolates with semi- 
axes (a,a,L), a <C L different power law behavior for 
parallel, perpendicular, and aligned configurations: 
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We have found no truly satisfactory approximation that 
would reproduce all intermediate limiting cases for el- 
lipsoids of arbitrary shape. Nevertheless we have made 
some progress compared to a simple 6-12 potential. Our 
somewhat naive strategy for developing an approxima- 
tion for the general case is (i) to use the Deryaguin ap- 
proximation in the short distance limit and (ii) to treat 
ellipsoids as spheres of equivalent volume on length scales 
which exceed the particle diameters. 



In the limit where large particles almost touch, the rel- 
evant distances become small compared to the local radii 
of curvature of the bodies. The short-distance expansion 
of Eq. (3) in powers of the local curvature radii is known 
as Deryaguin approximation [7,15]. In the most general 
case [16], each body has two different principal radii of 
curvature Ri and R' t at the point of closest approach. 
Furthermore, the principal axes of the two surfaces can 
be rotated by an angle 9 relative to each other. The re- 
sult of the Deryaguin approximation for Lennard- Jones 
interactions can be written in the form 
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where, in analogy to Eq. (5), we tentatively identify the 
orientation and relative position dependent part with a 
product of two terms to be specified below: 
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We note that Eq. (34) defines together with Eq. (28) 
and (29), a parameter-free 6-12 potential which correctly 
describes the position and depth of the energy minima 
for large colloidal particles. 

Our second step is the generalization of Eqs. (20) and 
(21) to the interaction between ellipsoids [17] 
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This ansatz reproduces at least the intermediate power 
laws for parallely and perpendicularly oriented thin pro- 
lates quite well: 
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In contrast, our ansatz overestimates the potential at in- 
termediate distances for needles aligned along their long 
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FIG. 3. Distance dependence of the potential energy for all possible pole contacts between an oblate and a prolate ellipsoid 
with semi-axes (1,6, 6)er and (2, 2, 9)a respectively: (•) Hamaker potential obtained via a Monte Carlo evaluation of the six 

dimensional integral Eq. (3); ( ) Deryaguin approximation Eqs. (33) and (40); (■ ■ ■) a GB 6-12 potential Eq. (5) adjusted 

to reproduce the position and depth of the minima in the Deryaguin approximation Eqs. (28) and (29); ( ) our proposal 

Eqs. (35) and (36). Two data sets were shifted along the j/-axis for clarity reasons. Note that none of the approximations 
contains freely adjustable parameters. 

The inset shows the ratio of the well depths of the Hamaker potential and of our potential as a function of the absolute well 
depth, (x) Results for pole contacts; (■) results for ellipsoids with randomly chosen relative orientations and positions whose 
distance is varied along the center-to-center line. A comparison of the two data sets allows an evaluation of the quality of our 
approximations for X12 and ni2 (see Section VI). 



axes. However, the following comparisons will show that 
the deviations are typically within a factor of two or three 
over the entire range of distances. 

V. POLE CONTACTS BETWEEN BIAXIAL 
ELLIPSOIDS 

The evaluation of Eq. (34) is straightforward for pole 
contacts between aligned ellipsoids with Ai = A2 so that 
9 = 0. For two biaxial ellipsoids which touch at their c 
poles, the principle curvature radii at the touch point are 

a 2 b 2 

Ri = and R' i = In this case Eq. (34) reduces to 



X12 »7ia = 




As in the case of interacting spheres, we can check the 
quality of our approximations by comparing to results 
obtained by, in the present case numerically, integrating 
Eq. (3) for pole contacts between oblate and prolate ellip- 
soids. Fig. 3 shows good agreement for the shape, depth, 
and position of the minima. In particular, the figure 
demonstrates that for small molecules Eqs. (35) and (36) 
provide a significant improvement over the Deryaguin ap- 
proximation (or a 6-12 potential taylored to reproduce 
the minima of the Deryaguin approximation). Figs, la 
and b show that we also reproduce the crossover to the 
asymptotic behavior at large distances reasonably well, 
at least relative to the parameter-free Gay-Berne poten- 
tial (Fig. lc). 
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FIG. 4. Deviations of approximate potentials from the 
true Hamaker potential as a function of the absolute value 
of the Hamaker potential: (a) using the correct distance 
of closest approach Eq. (7); (b) using the approximate GB 
distance function Eq. (8). We show results for contacts 
between an oblate and a prolate ellipsoid with semi-axes 
(a,b,c) — (1,6, 6)<7 and (a, b, c) = (2,2,9)a respectively and 
randomly chosen relative orientations and positions: (dark x ) 
our potential, (gray +) a Gay-Berne 6-12 potential adjusted 
to reproduce the energy minima within the Deryaguin aprox- 
imation. The solid lines show the corresponding results for 
pole contacts (see Fig. 1). 



VI. COMPUTABLE EXPRESSIONS FOR 
ARBITRARY CONTACT GEOMETRIES AND 
THE RELATION TO THE GAY-BERNE 
POTENTIAL 

In the next step, we take a closer look at the orientation 
dependent prefactor in Eq.(33). We are less interested in 
its exact calculation (i.e. the determination of the con- 
tact points, the local curvature radii and the angle be- 
tween the principle axis on the two surfaces) than in find- 
ing expressions which offer a good compromise between 
computability and correctness. Quite interestingly, our 
final expressions turn out to be closely related to those 
used in Gay-Berne potentials. 

We begin by providing definitions for X12 an d f]v2 
which, when multiplied with each other, reproduce 
Eq. (34) and whose structure resembles Eqs. (13) and 
(11): 
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Ki = \J Ri R[ is the Gaussian curvature at the touch 
point, £1(8) is a two dimensional rotation matrix. For 
two spherical surfaces with Ri = R[ and R2 = R' 2 the 
second term reduces to 7712 = 1. 

The GB definition Eqs. (13) - (15) of X12 agrees for 
pole contacts with Eq. (41) provided the parameter /j, is 
set to 
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and the three BFZ energy parameters Eq. (15) are iden- 
tified with the Gaussian curvatures at the three poles 
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The relation between the GB definition Eq. (11) of rj 12 
and Eq. (42) is less direct. The reason is that the cur- 
vature matrices in Eq. (44) characterize surfaces and, as 
a consequence, are two dimensional. Similarly, the an- 
gle 6 describes the relative rotation of the two surfaces 
around their common normal vector at the points of clos- 
est approach. In contrast, GB variables characterize the 
shape and orientation of three dimensional bodies. In the 
following, we will present a heuristic combination of GB 
variables which reproduces Eq. (42) for pole contacts. 

As a first step, we consider the definition of the Gaus- 
sian curvature in terms of the curvature matrix, Eq. (43). 
At the c-pole n ci can be written in the form 
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For arbitrary orientations it is tempting to replace by 
the projected diameter 



CTi(Ai,ri 2 ) = (fu AfS^Aj fi2 
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In this manner we arrive at the following approximative 
expression for 7712 : 

dettSJK + det^/,! 
(det[H 12 ]/(a 1+ a 2 )) 1/2 

H 12 (Ai,A 2) f) = —Aj S?Aa + -A^A 2 (51) 

<7l <7 2 

As a side result we note that the expressions for X12 and 
r/i2 can be further simplified for those contacts which 
dominate in the ordered phases of typical liquid crystals, 
i.e. contacts between similar poles of identical ellipsoids 
(01, h,ci) = (o 2 , b 2 , c 2 ) = (a, 6, c). 
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All in all, we are lead to expressions which are in remark- 
able agreement with those proposed by BFZ. In the case 
of xi 2 , we employ the same functional form together with 
particular choices for the four adjustable parameters /z, 
&ai, £bi, Sd- hi the case of 7712 , our Eqs. (50) and (51) re- 
spectively Eq. (53) resemble the corresponding Eqs. (11) 
and (10) so strongly, that it seems clear that our con- 
sideration eliminate with v = 1 the last remaining free 
parameter of the GB potential. 

We note that our proposed modifications leave with 
Xi 2 and ?7i 2 the most CPU time intensive part of the GB 
potential essentially unchanged. The small number of 
additional scalar operations necessary for the evaluation 
of Eqs. (35) and (36) hardly affects the performance of 
simulation codes. 



VII. NUMERICAL TEST OF THE 
APPROXIMATIONS FOR ARBITRARY 
RELATIVE POSITION AND ORIENTATION OF 
THE ELLIPSOIDS 

The most important question is, of course, how reliable 
the proposed approximations are. Similarly to Eqs. (35) 
and (36), the combination of Eqs. (13), (45), (47), (50) 
and (51) can only be considered as an educated guess 
for Eq. (34). The fact that we reproduce the results of 
the Deryaguin approximation for pole contacts inspires 
some confidence, but otherwise we have made substan- 
tial and uncontrolled approximations which need to be 
checked against the numerical evaluation of Eq. (3) for 
various relative positions and orientations. We represent 
the results by plotting the ratio of the approximative and 
the exact energy as a function of the exact energy. In 
this manner, results from a high-dimensional parameter 
space are (i) projected onto a single axis and (ii) sorted 
by importance. The inset in Fig. 3 shows a compari- 
son of the depths of the energy wells close to contact for 
fixed random orientations, while Fig. 4a deals with the 
attractive part of the interaction at arbitrary distances. 
Both figures also contain results for the pole contacts dis- 
cussed before to allow for an independent evaluation of 
the quality of the approximations for the distance and 
for the orientational part of the interaction potential. 

When judged against the corresponding results for the 
Gay-Berne potential, the agreement between our pro- 
posal and the numerical evaluation of Eq. (3) is excellent. 
In absolute terms, the deviations do not exceed a factor 
of two to three in cither direction. Quite interestingly, 
our approximations for X12 and ?7i 2 do not seem to be 
the source of large additional errors. Fig. 4b shows that 
the agreement is significantly reduced, if the approximate 



GB distance function Eq. (8) is used instead of the true 
distance of closest approach Eq. (7). 



VIII. SUMMARY 

We have presented an approximative interaction po- 
tential for soft ellipsoidal particles. Our potential uses 
(almost) the same variables as the Berardi, Fava and 
Zannoni [8] form of the Gay-Berne [3] potential for bi- 
axial ellipsoids, agrees significantly better with the nu- 
merically evaluated exact interaction potential, has no 
unphysical limits, and avoids the introduction of empir- 
ical adjustable parameters. The main modifications we 
propose are 

1 . to abandon the unphysical factorization of the ori- 
entation and distance dependent parts of the po- 
tential (Eq. 5) as well as the Lennard- Jones like 
form of the distance dependence itself (Eq. 6) and 
to replace them by Eqs. (35) and (36). 

2. to use the Gaussian curvatures at the ellipsoid poles 
Eqs. (45) - (47) in order to characterize the relative 
well depth for side-to-side, face-to-face and end-to- 
end interactions through the orientation and rela- 
tive position dependent factor xi 2 (Ai, A 2 , fi 2 ) de- 
fined in Eqs. (13) and (14). 

3. to replace the definition Eqs. (10) - (12) of the 
purely orientiation dependent factor r7i 2 (Ai,A 2 ) 
by Eqs. (49) - (51). 

4. to use the (effective) Hamaker constant A\ 2 [7] to 
set the energy scale. 

Our results for the attractive part of the soft potential are 
directly applicable to hard ellipsoids with van-der-Waals 
interactions. Furthermore, the proposed potential com- 
prises the interaction of point particles with ellipsoids as 
a well defined limit. This may be of interest for studies 
of wetting or polymer adsorption in colloidal dispersions 
which so far assume either a flat or a cylindrical geome- 
try [18]. 
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